//
// File:   main.cc
// Author: kajetan
//
// Created on 19 padziernik 2007, 19:43
//

#include <cstdlib>
#include <sstream>
#include <fstream>
#include <iostream>
//#include "model.h"
#include "algorytmyNumeryczne.h"
using namespace std;

string nazwaPliku = "y";
int dyn_disorder;

void drukuj(Model model, int st) {
    ostringstream oss;
    oss<<"output\\"<<nazwaPliku<<st;
    
    string c = oss.str();
    ofstream plik(c.c_str(), ios::out);
    cout<<c<<endl;
    
    double suma=0.;
    
    plik<<"# krok: "<<st<<endl<<"# czas: "<<st*DT<<" j.a.; "<<st*DT*2.41888e-17<<" s\n";
    for(int i=0; i<model.iloscBaz; ++i) {
        plik<<i<<" "<<model.getP(0,i)<<" "<<model.getP(1,i)<<endl;
        suma+=model.getP(0,i)+model.getP(1,i);
    }
    plik<<"# ---------------------\n# "<<suma<<endl;
    plik.close();
}

int main(int argc, char** argv) {
    
    string input, tmp;
    ostringstream oss;
    oss<<"output\\";
    int realizacji;
    
    input = (string)argv[1];
    realizacji = atoi(argv[2]);
    if(argc>4 && strcmp(argv[3], "-VACUUM")==0) {
        dyn_disorder = 1;
        tmp="VACUUM";
        oss<<"vacuum\\";
    }
    else if(argc>4 && strcmp(argv[3], "-STATIC")==0) {
        dyn_disorder = 0;
        tmp="NO";
        oss<<"static\\";
    }
    else if(argc>4 && strcmp(argv[3], "-WATER")==0) {
        dyn_disorder = 2;
        tmp="WATER";
        oss<<"water\\";
        cout<<"water";
    }
        
    for(int i=0; i<input.length(); ++i)
        oss<<(unsigned char)tolower((unsigned char)input[i]);
    ofstream plik(oss.str().c_str(), ios::out);
    oss<<".r";
    ofstream plik2(oss.str().c_str(), ios::out);
    oss<<".rk";
    ofstream plik3(oss.str().c_str(), ios::out);
        
    double tTot, sredni_czas = 0., odch = 0., tm, tLimit;
    int licz_realiz = 0;
    vector<double> times(realizacji, -1.);
    
    tLimit = (double)atoi(argv[4]);

    vector<double> sum(KROKOW, 0. );
    vector<double> part(KROKOW, 0. );    
    vector<int> freq(KROKOW, 0);
    
    plik<<" # multiplier="<<multiplier<<"\n # num. | t | surv\n";
    for(int r=0; r<realizacji; ++r) {
        Model model(input);
        model.inicjuj();
        
        if(argc>5)
            model.GAMMA = (double)atoi(argv[5])/1000000.;
        else 
            model.GAMMA = 1.e-4;
        
        tm = -1.;
        tTot = 0.;
        freq.clear();
        part.clear(); 
        
        int i, s;
        double suma;
        
        for(i=1; i<=(KROKOW/KROKOW_NA_WYDRUK); ++i) {
            algNum::iteruj(KROKOW_NA_WYDRUK, model, dyn_disorder, tTot);
           
            suma = 0.;
            plik3<<tTot<<" ";
            for(int j=0; j<model.iloscBaz; ++j) {
                plik3<<model.getP(0,j)<<" ";
            }
            for(int j=0; j<model.iloscBaz-3; ++j)  {
                suma += model.getP(0,j);
                suma += model.getP(1,j);
            }
            plik3<<endl;
                
            if(tm < 0. && suma * Em < 1.) {
                tm = tTot;
                sredni_czas += tTot;
                ++licz_realiz;
                times[r] = tTot;
            }

            if(tTot > tLimit)    break;
                       
            s = (int)(tTot / DT);
            part[s] += suma;
            ++freq[s];
        }
        suma=0.;
        if(model.getP(0,0)>1.) {
            plik<<r<<" "<<tTot<<" -1. : "<<model.getP(0,0)<<endl;
            for(int i=0; i<model.iloscBaz; ++i) {
                plik<<"\t"<<i<<" "<<model.getP(0,i)<<" "<<model.getP(1,i)<<endl;
            }
            break;
        } else {
            for(int j=0; j<model.iloscBaz-3; ++j)
                suma+=model.getP(0,j)+model.getP(1,j);
            plik<<r<<" "<<tTot<<" "<<tm<<" "<<suma<<endl;
        }
        
        for(int j=0; j<KROKOW; ++j) {
            sum[j] += (!freq[j]) ? 0. : part[j]/freq[j];
        }
        cout<<model.iloscBaz<<endl;
    }

    sredni_czas /= licz_realiz;
    for(int r=0; r<realizacji; ++r) {
        if(times[r]>=0.) {
            odch += (sredni_czas-times[r])*(sredni_czas-times[r]);
        }
    }
    odch /= licz_realiz;
    
    plik<<"# average: "<<sredni_czas <<" variation: "<<sqrt(odch)<<endl;
    plik<<"# complete: "<< licz_realiz <<" incomplete: "<<(realizacji - licz_realiz)<<endl;
    
    for(int j=1; j<KROKOW && (!j || sum[j]>1.e-2); ++j) {
        plik2 << j*DT << " " << sum[j] / realizacji << endl;
    }
    
    plik.close();
    plik2.close();
    plik3.close();
        
    return (EXIT_SUCCESS);
}
